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We present a theoretical analysis of a coupled, two-state Bose-Einstein condensate with non- 
equal scattering lengths, and show that dynamical instabilities can be excited. We demonstrate that 
these instabilities are exponentially amplified resulting in highly-directional, oppositely-propagating, 
coherent matter beams at specific momenta. To accomplish this we prove that the mean field of 
our system is periodic, and extend the standard Bogoliubov approach to consider a time-dependent, 
but cyclic, background. This allows us to use Floquet's theorem to gain analytic insight into such 
systems, rather than employing the usual Bogoliubov-de Gennes approach, which is usually limited 
to numerical solutions. We apply our theory to the metastable Helium atom laser experiment of Dall 
et al. [Phys. Rev. A 79, 01 1601 (R) (2009)] and show it explains the anomalous beam profiles they 
observed. Finally we demonstrate the paired particle beams will be EPR-entangled on formation. 
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I. INTRODUCTION 

The experimental realization of Bose-Einstein conden- 
sates (BECs) [TJ [2] has allowed the creation of macro- 
scopic quantum systems that are sensitive, isolated from 
the environment, and highly controllable, providing an 
excellent system in which to investigate many body quan- 
tum mechanics [3HB] , as well as leading to the creation of 
the atom laser, a matter-wave analog of the optical laser 

Atom lasers are usually generated by creating a BEC 
in which the atoms are trapped by magnetic or optical 
fields, and then changing the state of a subset of the 
atoms so that they are no longer trapped, and are free to 
escape. This can be done via RF spin flipping [7], or Ra- 
man outcoupling [HI QUI , respectively resulting in a beam 
that falls under gravity, or a beam that can have some 
directionality due to a momentum transfer correspond- 
ing to a few cms -1 . Such coherent, directional matter 
beams offer great promise not only as a sensitive probe of 
the source condensate, but also offer potential in terms 
of atom interferometry and precision sensors. Further- 
more, both BECs and atom lasers also offer the possibil- 
ity of creating non-classical states that exhibit properties 
such as number and quadrature squeezing [TTrHS] , which 
among other things have the potential to enhance the 
sensitivity of atom interferometers |16j , allowing the mea- 
surement of electric, magnetic and gravitational fields, as 
well as rotations, with great precision [TT1EU] . 

The creation of coherent matter beams can also be 
achieved by utilizing the inherent properties of the BEC 
itself. This can be more useful than a conventional atom 
laser as not only can such beams carry information about 
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the interaction properties, dynamics and quantum state 
of the BEC, but depending on their method of genera- 
tion they can possess highly interesting quantum states 
of their own. The dissociation of a molecular BEC has 
been shown, for example, to lead to the creation of pairs 
of particles with oppositely-directed momenta that are 
expected to exhibit EPR correlations [27], as can differ- 
ent forms of atomic four- wave mixing [28 31j . One can 
also utilize colliding condensates resulting in a spheri- 
cal shell of scattered atoms where opposite points on the 
sphere have correlations that depend on the the intera- 
tion properties of the atoms within the condensate. Most 
of these schemes, however, result in atomic fields that are 
not in easily accessible modes. They rely, for example, on 
collisions or dissociation resulting in scattering over the 
full sphere, beams that are co-propagating rather than 
spatially distinct, or that have their interesting quantum 
properties such as entanglement spread over a large num- 
ber of resonant modes. 

In this paper we present an analysis of a new method of 
generating correlated atomic fields, which has the advan- 
tages of being experimentally realistic and producing spa- 
tially distinct, coherent, highly-directed atomic beams at 
well-defined momenta. The beams are produced due to 
dynamic instabilities within the atomic system itself, and 
as such do not rely on molecular dissociation or colliding 
condensates, or the existence of engineered quantum op- 
tical fields to use as a template for the desired quantum 
statistics to be transferred to the atomic fields. 

Our scheme is motivated by the experiment of Dall et 
al. [32] . who exploited the combination of a metastable 
helium BEC and a double stacked multichannel plate to 
image the transverse (i.e. orthogonal to the propaga- 
tion direction) beam profile of an atom laser. At high 
outcoupling powers they observed an anomalous particle 
flux that was ejected along the long axis of the BEC, 
and was well separated from the bulk atom laser profile. 
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In this paper we develop a general theoretical model of 
a class of systems which includes that experiment, and 
demonstrate that pairs of correlated beams at specific 
quantized momenta can be produced. 

In brief, our scheme requires a two-state BEC that has 
a high aspect ratio, with the two states being coupled in 
some way. A further requirement is that the s-wave scat- 
tering length of the atoms within each state and between 
the states must not all be equal. One experimental real- 
ization of such a system would be a magnetically trapped 
condensate with two equally-trapped Zeeman states, for 
example the F = 2, mp = 1 and F = 1, mp = —1 
levels of 87 Rb. Another would be an optically trapped 
condensate using the quadratic Zeeman shift to lift the 
degeneracy in the Zeeman level splittings to create an 
effective two-level system. In both the optical and mag- 
netically trapped cases the two levels could be coupled 
by two optical fields utilizing a Raman transition. Two- 
state Raman coupling has been achieved experimentally 
in 87 Rb [23. 

The existence of the coupling between the two levels, 
as well as the non-equal scattering lengths, results in 
dynamic instabilities in the condensate which can grow 
exponentially at specific well-defined momenta. As we 
will demonstrate, the equations governing the unstable 
system can be mapped to those describing optical para- 
metric down-conversion, and have the same result — the 
generation of pairs of entangled particles with oppositely 
directed momenta. Only the unstable momentum modes 
which have de Broglie wavelengths that fit within the 
condensate are amplified, which results in the pairs be- 
ing emitted along the long (weakly trapped) axis of the 
condensate. If the aspect ratio of the condensate is suf- 
ficiently large, a pair of highly directed beams will be 
generated. 

As part of our analysis, we will develop techniques for 
the application of Bogoliubov theory to perturbations 
about a time-dependent mean field background, rather 
than a static one. A Bogoliubov approach is necessary 
because the creation of the entangled pairs of particles 
is spontaneously seeded, meaning that analysis beyond 
the mean field is required. In its standard form Bogoli- 
ubov theory assumes a time-independent mean field. In 
our case, however, the effect depends on the existence 
of a time- dependent mean field background, which sug- 
gests a Bogoliubov-de Gennes approach should be used. 
The usual drawback to such an approach is that the 
Bogoliubov-de Gennes equations are generally analyti- 
cally intractable and are therefore solved numerically, 
which in many cases precludes deep insight into the prob- 
lem. In this paper, however, we prove that the mean field 
background in our system, while time-dependent, is also 
periodic. This enables us to employ Floquet's theorem 
to gain a semi-analytic solution to the problem of deter- 
mining the excitation spectrum and instabilities of the 
system. 

The outline of this paper is follows: We begin by giv- 
ing a brief overview of Bogoliubov theory in Section [TTl 



and demonstrate that the basic approach is still sound 
even in the presence of a time-dependent background. In 
Section |III| we solve for the mean field dynamics of our 
system, and map them to modified optical Bloch equa- 
tions with a nonlinear term. This enables us to prove that 
the mean fields, while time-dependent, are also periodic, 
regardless of the initial state of the system. The equa- 
tions of motion for the fluctuations beyond the mean field 



are derived in Section IV and we show that in the spe- 



cial case of equal nonlinearities for the two atomic states 
they reduce to the standard Bogoliubov excitation spec- 
trum. In Section [V] we take advantage of the temporal 
periodicity of the mean field to employ Floquet's theo- 
rem and show that the Floquet exponents can indicate 
dynamical instabilities in the system. These dynamical 
instabilities lead to exponential growth of excitations at 
particular momenta. We also check the accuracy of the 
theory by performing truncated Wigner simulations of 
the system to ensure the excitations occur where our Flo- 
quet theory predicts, and derive the equations of motion 
for the annihilation and creation operators governing the 
excitations. As a further test of the validity of our the- 
ory we apply it to the experiment of Dall et al. 32 in 
Section [Vl] and demonstrate that the anomalous particle 
production they measured is correctly described by our 
theory. In Section |VII| we note that the equations of mo- 
tion for the excitations are identical to those governing 
optical parametric down conversion, and thus should ex- 
hibit EPR entanglement, at least on formation. Finally, 
in Section [VIII[ we discuss the optimal experimental sys- 
tem to generate the dynamical instabilities our theory 
predicts, and consider possible difficulties, and conclude 
in Section Hxl 



II. OVERVIEW OF BOGOLIUBOV THEORY 

As described in the Introduction, the paired matter 
beams generated in our scheme arise from a dynamical 
instability in the condensate, and in order to determine 
which modes are dynamically unstable, the excitation 
spectrum of the condensate must be obtained. To that 
end, we begin with an overview of the theory of conden- 
sate excitations beyond the mean field approximation, 
and then use these results to calculate the stability of 
the condensate and its excitation spectrum in Section 
|IV| More comprehensive treatments of the Bogoliubov 
theory can be found in [35] and in a number of review 
articles (see, for example, [36U38] ). 

The problem is to determine the response of the con- 
densate to small fluctuations about the mean field. Typ- 
ically the condensate is stable to such fluctuations, and 
their energy spectrum determines the phase and group 
velocities of the excitations. In the case that the con- 
densate is dynamically unstable some modes will un- 
dergo exponential growth, which corresponds to the gen- 
eralised energy spectrum containing non-zero imaginary 
components. As a concrete example of the techniques 
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we will employ, consider a single-component condensate 
described by the Hamiltonian 



H 
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where an arbitrary energy offset fj, has been included. 
This term is introduced for calculational reasons and has 
no physical influence on the Hamiltonian as it cannot af- 
fect any expectation values. This can be seen by noting 
that although it contributes a different energy offset for 
states with different total number, these states are un- 
coupled due to the Hamiltonian conserving particle 
number. For the case of calculating the excitation spec- 
trum of a ground state, the arbitrary energy offset \i will 
be identified as the condensate chemical potential. 

To find the excitation spectrum of Eq. ([IJ we use the 
Gross-Pitaevskii equation to describe the mean field of 
the condensate, and then consider quantum-mechanical 
fluctuations about this mean field. To this end we define 
the deviation operator S^f = $ — where ^ = and 
S^f is considered to be a small quantity in the sense that 
its quadratic moments are small compared to |\&| (its 
first-order moments are identically zero). 

The concept of S^f being a "small quantity" is not en- 
tirely trivial. For it to be considered in any sense small, 
the mean field (Jff\ must be non-zero. However the state 
of condensates with a large number of atoms is well ap- 
proximated by either a state with well-defined total num- 
ber or as a statistical mixture of coherent states with 
random phases [3S]. In both of these cases the mean 
field (^) is zero, as the system has no well-defined global 
phase. As any physical expectation value is indepen- 
dent of the choice of global phase, however, any anal- 
ysis can be performed for a coherent state with a spe- 
cific global phase, and the results then averaged over the 
global phase. Because all observables will be made up 
of equal numbers of annihilation and creation operators 
(since atom number is conserved), the choice of global 
phase cannot affect the final result. 

An alternative justification is to note that spontaneous 
symmetry breaking results in the selection of a specific 
phase, which is equivalent to assuming that the state of 
the condensate is unchanged on the addition of a par- 
ticle, that is the states \N) and \N + 1) are physically 
identical [40l [41] . This is a robust assumption for BECs 
with particle number much greater than one, which is the 
case for most bulk condensate experiments. 

Making the replacement ^ — > ty + Sty, the Hamiltonian 
Eq. ill can be expanded in powers of <5^, 



H = H + H 1 + H 2 + H 3 + #4, 



(2) 



where H n contains terms of order (<5W)™. The excitation 
spectrum of the Hamiltonian H is then approximately 
given by the eigenvalue spectrum of the lowest order non- 
trivial term in Eq. ([2]). 



The zeroth order term in Eq. |2]) 
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is simply a constant and represents the total energy of 
the unexcited mean field. The first order term is of the 
form 



,/ x j m?pj 6%, (4) 



where the mean field \& evolves as 
dt ~ { 2M 



gib / ti 2 X7 2 
ih~= l-^r^+V(x) + U\y -/*)*. (5) 



Although the first order term H\ is non-zero, it does 
not affect the evolution of the deviation operator. To see 
this we note that as *f? = $ + 5$ we have 

/9* 3 - - - ~ 

ih—+ih-5y=[*,H] + [5*,H] (6) 

and because is a complex-valued function it commutes 
with all operators, giving 

ih— 5V = \8%H] - th- 
at 



dt 



= [<J* , + [8%, H 2 + H 3 + H A ] - ih 

d$> r ~ ~ - - i 9* 
ih— + [<J* , H 2 + H 3 + H 4 ] - ih- 



d^ 
dt 



dt 



dt 



(7) 



Since Hi does not occur in Eq. <[7j) , it does not affect the 
evolution of the deviation operators and will not con- 
tribute to the excitation spectrum. We note that typi- 
cal treatments of the Bogoliubov theory consider the re- 
stricted case of a static condensate density and choose 
/i as the chemical potential such that d^ /dt = 0, and 
hence Hi = 0. As shown by Eq. this choice of /i 
is unnecessary as Hi does not influence the evolution of 
£>]/, independent of the choice of /i and even in the gen- 
eral case of a non-stationary mean field. This is relevant 
as we will be considering the extension of this theory to 
a condensate with a time-dependent density background. 
The first physically important term in Eq. |2]) is 



Ho = 



- / h 2 X7 2 9 
dx8& + V(x) + 2U\t>\ - a* ) 5* 

hj J dx (V<S*W + (**) 2 <5*<5* 



(8) 



As the deviation operator 8^ is small compared to the 
mean field \P, the higher-order contributions H 3 and H^ 
to the total Hamiltonian can then be neglected compared 
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to H 2 . The excitation spectrum of the condensate about 
the mean field \& is then approximately given by the en- 
ergy spectrum of H 2 . 

To avoid solving the infinite dimensional eigenvalue 
problem H 2 \^) = E\ty) for the condensate excitation 
spectrum directly, it is desirable to apply a linear trans- 
formation to H 2 that will diagonalize it in the form 



h 2 = J2^1K 



(9) 



for some boson annihilation operators Aj and real fre- 
quencies LUi. In this form, the Hamiltonian can be simply 
interpreted as representing a set of modes with energies 
hiOi, which is the condensate excitation spectrum. In 
general, it is not possible to transform H 2 into the form 
Eq. Q if the Hamiltonian possesses any instabilities [12] , 
however one frequently considers the excitation spectrum 
of the ground state which is stable by definition and in 
this case such a transformation is always possible. 

In the general case, we look for the operators Aj satis- 
fying 



d ~ 

ih— Ai = [Ai,H 2 ] = -hwiAi 



(10) 



where Wj is real if and only if Aj is a boson annihilation 
operator |42| . In the case that uji is complex, boson anni- 
hilation operators can be constructed from the Aj as dis- 
cussed in Appendix|Aj Hence huii can be considered to be 
a generalised energy spectrum of the condensate where 
non-zero imaginary components correspond to dynami- 
cal instabilities of the condensate. Note that although 
the eigenvalues of H 2 must be real as it is Hermitian, the 
eigenvalues of Eq. ( 10 ) need not be real. For example, the 



Hamiltonian for degenerate parametric down-conversion 
H = hx (aa + a 1 a^) is Hermitian but the correspond- 
ing eigenvalues of Eq. ( 10 ) are pure imaginary. In this 



case the occupation of the mode a undergoes exponential 
growth. The case of complex eigenvalues uji is discussed 
further in Sections |IV| and |VJ and Appendix |Aj 

Equation ( 10 ) is most easily solved by expanding the Aj 



in a complete, linearly independent basis {Tj} such that 



A, 



cjY where Cj is a complex vector, cj denotes its 



conjugate-transpose, and Y is the column vector formed 



by the complete basis {Tj}. For the case of Eq. |8|), an 
appropriate basis is {Tj} = {S^>, S^ 1 ^}. As the Hamilto- 
nian H 2 is quadratic, its commutator with every opera- 
tor Tj will be linear in the operators {Tj}. Defining the 
complex matrix % to represent this relationship 



£> Jfc T fc = [Tj,^] 



(11) 



permits Eq. ( 10 ) to be recast as an eigenvalue problem 
in "H, 



[cl±,H 2 ] = <$UY = -fUJiclt, 



-hioic] 



(12) 
(13) 



where the last line follows as the components of Y 
are linearly independent. If the mean field VP is time- 
independent, then the matrix % will also be time- 



independent and Eq. ( 13 ) represents an eigenvalue prob- 



lem for the left eigenvectors cj and eigenvalues — HtJi of 
the matrix H. If the mean field simply evolves due to 
a global phase rotation, this can be cancelled by appro- 
priate choice of the arbitrary energy offset fi making % 
time-independent. In the case of the condensate ground 
state, that offset will be the chemical potential of the 
condensate. The eigenvalues {hu>i} then represent the 
generalised excitation spectrum of the condensate about 
the mean field, which was to be determined. 

It is important to note that for the eigenvalues of H. 
to determine the solution to Eq. ( 10 ) the matrix % must 
be constant. In the next section these techniques will be 
generalised to the case of a periodic mean field in which 
the time-dependence of the matrix % cannot be removed 
by any analytic transformation. 

In the case of a homogeneous condensate (V(x) = 0) 
the matrix H can be diagonalized analytically to give the 
condensate excitation spectrum as 



hu>(k) = y/e(k) (e(k) + 2nU), 



(14) 



where k is the wave vector, e(k) = h 2 k 2 /2M is the free- 

1 1 2 

particle energy spectrum and n = \fM is the condensate 
density. Equation |l4| is known as the Bogoliubov exci- 
tation spectrum |43j. 

The interested reader is referred to the review paper 
by Ozeri for further details about Bogoliubov theory in 
Bose-Einstein condensates [37] . 



III. MEAN FIELD DYNAMICS AND 
EXISTENCE OF PERIODIC SOLUTIONS 

As described in the Introduction, our scheme consists 
of a two- level condensate with the two levels |0) and |1) 
coupled in some way, leading to Rabi oscillations between 
the levels. 

Under this approximation, the second-quantized 
Hamiltonian for the system is given by 

i=l,0 J 



2 E u *i I d**I*}*i* 4 

a a in " 



(15) 



ij=l,0 



+ nn I dx (* t 1 # + *J ) * 1 



where Uij — Ai:h 2 aij/M is the nonlinear interaction 
strength, ajj is the s-wave scattering length between in- 
ternal states \i) and f2 is the Rabi frequency which 
is taken to be real, and fi is an energy offset which has 
been included to cancel the global phase rotation which 
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would otherwise be present. The equations of motion 
corresponding to this Hamiltonian are 

+/ifi#o-A»*i, (16) 



ih— *o = -2^ V 2 # +C/(* t 1 * 1 + K^* ) * 



dt 



-Ml^x - /i* , (17) 



where U = U\\ = Uio and we have defined the nonlin- 
earity mismatch parameter 



K = Uqo/Uh. 



(18) 



As described in Section|TTJ the excitation spectrum of a 
condensate can be obtained by approximating each field 
operator as a complex number (the mean field) plus a 
small fluctuation term, and then either diagonalizing the 
Hamiltonian [43, 44 a or diagonalizing the linearized equa- 
tions of motion for the fluctuations themselves. We will 
take the latter approach in this paper, but with the dif- 
ference that the mean field about which the linearization 
procedure will take place is itself time- dependent. 

The method for diagonalizing the evolution equations 
of the linearized fluctuations to obtain the excitation 
spectrum about the ground state is the same method 
used to determine the stability of fixed points of systems 
of nonlinear ordinary differential equations. As men- 
tioned in the previous section, this method relies criti- 
cally on the fact that it is a stationary solution about 
which the equations are linearized. Floquet's Theorem 
[15] allows the stability of periodic solutions to be con- 
sidered, and it is this theorem that will be used to de- 
termine the stability of the condensate about these mean 
field dynamics. We will now show that the mean field 



evolution of Eq. (17 1 is periodic 



Within the local density approximation [JS] [17], we 
can assume that the mean field remains homogeneous; 
only the excitations will have spatial dependence. The 
equations of motion for the mean field then reduce to the 
following ordinary differential equations, 



= U (|*i| 2 + |* | 2 ) * i + KM - /i*!, (19a) 
ihj^o = U (|*i| 2 + k|* | 2 ) *o + fiOtfi - m*o- (19b) 



Although solving Eq. ( 19 ) for k ^ 1 is intractable ana- 



lytically, it can at least be shown that the solutions are 
periodic up to a global phase rotation, and exactly peri- 
odic for appropriate choice of the arbitrary energy offset 
fj,. The simplest way to do this is to recognise that these 
equations are modified optical Bloch equations contain- 
ing a nonlinear term but with no damping. Defining 
*i = Ciy/n where n = |*x| 2 + |* | 2 is the total den- 
sity, the equations of motion for the density matrix terms 



Pio = ^Cq and w = p n - p 0Q = | c x | 2 - |c | 2 are 



d pio = -if (1 - w)p w + iQw, 



dt 



dt 



-4fflmf>io), 



(20a) 
(20b) 



where g = nU(l— k)/H. As the evolution is purely Hamil- 
tonian, these equations conserve the expectation value of 
the Hamiltonian. In particular, as the number of atoms 
is also conserved, the mean energy per particle given by 



E 



(Ml 
(N) 



-M 1 



+ 2Mme(pi ) (21) 



is conserved. The solutions to Eq. (20 1 are plotted in 
Figure [T] 

As the evolution is purely Hamiltonian, the state can 
be described by a point on the surface of the Bloch sphere 
(see Figure [l]). The state is however not completely 
free to move on this sphere as the energy per particle 
E must be conserved by its motion. Equation ( plj ) is 
a holonomic constraint, which together with the identity 
w 2 + 4|/9io| 2 = 1 reduce the number of degrees of freedom 
of the solution from three (w, 2 Re(pio), 2 Im(pio)) to one, 
restricting the system's motion to a one-dimensional sub- 
set of the surface of the Bloch sphere (lines of constant 
color in Figure [lj. 

There exists an extension to the Poincare-Bendixson 
theorem which states that for certain manifolds, includ- 
ing the sphere S 2 , all possible paths traced out by a C 2 
action (which includes paths in phase space from a suffi- 
ciently smooth Hamiltonian) must approach or be peri- 
odic cycles, go to fixed points, or be a curve that covers 
the entire surface [IS] . On a sphere the last possibility is 
clearly impossible due to the hairy ball theorem demon- 
strating the existence of at least one fixed point [49] . 

In Appendix[B]we show that the fixed points in the cur- 
rent system are unreachable, meaning they either cannot 
be reached in normal evolution, or are a trivial case of 
periodicity where nothing changes. The final possibility 
consistent with the theorem proved in [15] is that some 
trajectories may approach a limit cycle. In Appendix [B] 
we also show that limit cycles cannot exist in this sys- 
tem. Consequently, all trajectories in this system must 
be periodic. 

Although it follows from the periodicity of the evolu- 
tion of the state on the Bloch sphere that all physical 
expectation values are periodic (as they must be inde- 
pendent of the global phase), the global phase rotation 
must be cancelled for the wave functions themselves to 
be periodic. This can be achieved by appropriate choice 
of the energy offset p. It is the periodicity of the or- 
der parameters that will enable the stability of the 
condensate to excitations to be determined. 
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FIG. 1: (color online) Bloch sphere representation of the evo- 
lution described by Eq. (201. The upper figures (a) and (b) 



represent the case of the usual optical Bloch equations with no 
damping (g/Q = 0), middle figures (c) and (d) illustrate the 
effect of the nonlinear term on the evolution with g/Q — 8, 
and the lower figures (e) and (f) illustrate the limit where the 
nonlinear term dominates (g/Q — > oo). The left figures (a), 
(c), (e) illustrate the Bloch sphere colored according to the en- 
ergy per particle E (see Eq. (21 1). The system is constrained 
to move on lines of constant color. Bands have been removed 
from these spheres for illustration purposes only. The right 
figures (b), (d), (f) are contour plots of the energy per particle 
E over the surface of the Bloch sphere. 



IV. EXCITATION DYNAMICS 

The evolution of small perturbations about the mean 
field dynamics of a condensate define both the excitation 
spectrum of the condensate and its stability to perturba- 
tions. 

To determine the evolution of these excitations the 
mean field dynamics must be separated from that of the 
excitations. As in Section III] we proceed by defining the 



deviation operators S^i = ^ — y&i) and treat as 
a small quantity. In this case, however, the (JjM = &i 



are themselves time dependent, obeying Eqs. ( 19 1 for the 
mean field. 

The equations of motion for the deviation operators are 
obtained by replacing the field operators ^ in the oper- 
ator evolution equations (17) with ^ +5^i and keeping 



only terms up to first order in the deviation operators. 
Applying this procedure gives 



at 



=U 



(22a) 



2M 



V ifi + hfLS^o - A*<5*i, 



at 



--U 



(2 K |* | 2 + |* 1 | 2 ) <y* + /t*g^t 



+*1 ^qJ^J + ^^o^^/i 

h 2 



(22b) 



2M 



V"W + - [iS^o- 



Having assumed the mean field (but not the fluctua- 
tions) to be homogeneous, the evolution equations are 
spatially translation-invariant and take their simplest 
form in a Fourier basis. Performing the Fourier trans- 
form of Eqs. (22) yields 



iftj^*i(k) 



U 



(2|* 1 | 2 + |* | 2 )^i(k) 



+^iS<b[(-k) + *!* ^(-k) 
h 2 k 2 



Vt, 



+*5*i<5*o(k) 



-<J*i(k) 



ih-s^ (k) 



= u 



2M 

+m8V (k) - fi6$i (k) , (23a) 

0| 2 + |* 1 | 2 )^ (k) 



+/c*g«*S(-k) + *i*o«*I(-k) 
S^o (k) 



Vt/ 



-*i*o^i(k)j 2V/ 
-MM#i(k) -//5* (k) 



(23b) 



In this form, it is clear that the Fourier modes 
are almost completely decoupled from each other. 
Each deviation operator 5^i(k) is only coupled to 

(s4fj(k), f>*j(-k) j, with each <S*J(— k) also only cou- 
pled to this same set. This can be exploited to write 



Eqs. (231 in matrix form as 



ift^*(k) = w(k)f (k) 



(24) 



where 



t(k) = (^(k) S<t>\(-k) S^ (k) S¥ (~k)) , (25) 
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and gi = £/(2|*i| 2 



<Zo 



C/(2ac|* | 2 



l*i| 2 ), 



ii SJ = ?7*?*j, «ij = f7#i% and e(k) = h 2 k 2 /2M 

Note that the matrix H(k) is not the Hamiltonian, but 
is related to it by Eq. (111. As a consequence, although 



it will be shown later that in some circumstances H(k) 
contains complex eigenvalues and is hence not Hermitian, 
this in no way conflicts with the requirement that the 
Hamiltonian H must be Hermitian and only have real 
eigenvalues. 

If the coefficients of the matrix 'H(k) were not time- 
dependent, the excitation spectrum of the condensate 
could simply be obtained from the eigenvalues of H(k). 
Non-zero imaginary components for these eigenvalues 
would indicate the corresponding mode to be unstable. 
This is true independent of the sign of the imaginary 
component as the eigenvalues come in pairs with oppo- 
site imaginary components. A full discussion of this issue 
is given in Appendix |A"| 

Before considering the general case of k ^ 1, it is 
instructive to consider the limit in which all scattering 
lengths are equal (k — 1) and recover some familiar re- 
sults. If we assume k = 1, the nonlinear term in Eq. ( 19 ) 



only contributes to a rotation of the global phase of the 
spinor condensate. In this case the dynamics can be 
solved analytically and familiar excitation spectra recov- 
ered. We find that the general solution is 

#i(t) = cos(fii)$+ + sin(fi*)$_, (27a) 
#o(*) = -isin(fit)$ + + £cos(fit)$_, (27b) 

for some complex constants $± , and where the chemical 
potential /i = nil has cancelled the global phase rotation. 
This solution can be viewed as a linear basis transforma- 
tion from to <I>±, which are the eigenvectors of the 
Rabi coupling term in Eq. ( |15[ ). Performing this change 
of basis on the original Hamiltonian given by Eq. ( 15 ) 



yields a Hamiltonian of the same form, but without the 
Rabi coupling term. The equations of motion for the de- 
viation operators 6&± ther efore give a matrix of precisely 
the same form as Eq. ( 26 ) but in terms of <I>± and SQ± 
instead of the VP^ and 5^>i, and with f2 replaced by 0. 
This new matrix "H'(k) is time-independent and can be 
diagonalized to give the eigenvalues 

fedt(k) = y/e{k.) (e(k) + 2nC/), (28a) 
fe^(k) = e(k), (28b) 



where e(k) = h 2 k 2 /2M , n 
maining two eigenvalues are the negatives of Eq. ( 28 1 



hu>f(k) is the usual Bogoliubov spectrum |43j correspond- 
ing to excitations in the total condensate density. The 
eigenvalue So;j.(k) is the free particle spectrum; this ex- 
citation only changes the relative densities of the two 
states without affecting the total density, hence not af- 
fecting the nonlinear term in the Hamiltonian given by 
Eq. g5j. 

As explained in Section [TTJ the Hamiltonian for the 
condensate excitations that corresponds to the eigenval- 
ues given by Eq. ( 28 ) is 



H = 



(29) 



where the A^ ^(k) obey boson commutation relations 
and are the corresponding normalised eigenvectors to the 
eigenvalues in Eq. ( 28 1 . The normalised eigenvectors for 



the negatives of those eigenvalues are the A| 4 (-k). 



V. INSTABILITIES AND THE EXCITATION 
SPECTRUM 

Having considered the limit of equal scattering lengths, 
it now remains to determine the energy spectrum and 
condensate stability in the general case of k ^ 1. 

In the general case, the excitation spectrum cannot 
be obtained from the eigenvalues of the matrix H(k) 
in Eq. (26) as the matrix's entries are themselves time- 
dependent. However, due to the periodicity of the mean 
field wave functions demonstrated in Section IIIIl the en- 



tries of "H(k) are themselves periodic, which enables Flo- 
quet's theorem to be applied. 

Floquet's theorem proves that the matrix solution to 
the initial-value problem 



-U(t) = A(t)Tl(t), 

n(o) = i, 



(30a) 
(30b) 



where I is the n x n identity matrix and A(t) a periodic 
n x n matrix with period T can be written in the form 



II(i) = P(t)exp(-iQt), 



(31) 



where Q is some constant matrix, P(t) is a matrix of 
periodic functions with period T, and P(0) = I. 

We note there are differing definitions of the matrix 
Q. While it is usual in quantum mechanics literature 
to define Eq. (31) with the — i in the exponent, 



in mathematics literature the —i is omitted (see, for ex- 
ample ma Eg). 

The matrix solution H(t) is the general solution to the 
related linear system 



dt 



x(t) = A(t)x(t), 



(32) 



for any initial condition x(0) where x(t) is a vector. Ev- 
ery solution x(t) to this problem can be written in terms 
of the matrix Tl(t) using 



x(t) = n{t)x(o), 



(33) 



as is easily verified. The matrix solution Tl(t) thus 
completely determines the behaviour of all solutions to 
Eq. ((321), or, equivalently, Eq. E4h. 



The eigenvalues of Q are known as Floquet exponents 
(or characteristic exponents) and determine the long- 
term growth or decay of the solutions to Eq. ( 32 ) . These 



eigenvalues can be obtained from the monodromy matrix, 



M = n(T) = cxp(-iQT), 



(34) 



as P{T) = P(0) = I. The existence and uniqueness 
of the solution to Eq. (30 1 guarantees that U(t) and 
hence M. will be invertible. The Floquet exponents £j 
can therefore be obtained from the eigenvalues Xi of the 
monodromy matrix using = exp(— z£,T). It is the 
Floquet exponents of the matrix Ti(k) that we wish to 
calculate in order to determine the stability of the con- 
densate to excitations, as these exponents will determine 
whether the long term behaviour of the system exhibits 
exponential growth or decay. 

Determining the Floquet exponents of the system via 
Eq. ( 24 1 requires knowledge of its period T, and hence 



the period of the mean field dynamics given by Eq. ( 19 1 



Although this period cannot generally be determined an- 
alytically, it can be found numerically. 

It was shown in Section |Hl| that up to a global phase 
rotation f(T) = e i27rA " T /(0), the mean fields are 
periodic. The mean fields ^j(t) can be therefore written 
in the form 



^j(t) = aj iTl exp [i2n (nv + Av) t] , 



(35) 



for some complex constants <x/,„, fundamental fre- 
quency vq = T^ 1 , and frequency offset Av. In this 
form, the \& j(t) are not exactly periodic as *f?j(T) — 
exp(iAi/r)^ j (0) , but this frequency offset can be can- 
celled by an appropriate choice of the energy offset /x = 
-2-kKAv in Eq. (pi 



The period T and frequency offset Av in Eq. (35) can 



be determined from the Fourier transform of \& j (t), which 
will have sharp peaks at the frequencies jivq+Av. Choos- 
ing the energy offset fi = —2-kTiAv, the frequency offset 
in Eq. ( 35 1 can be cancelled making the tyj (t) with this 
energy offset exactly periodic with period T. 




Frequency (kHz) 



FIG. 2: (color online) Temporal Fourier transform of the cal- 
culated mean field evolution defined by Eq. ( |19[ l for the centre 
of a two-state metastable helium condensate N — 2 x 10 6 , 
uj r = 2tt x 1020Hz, uj z = 2tt x 55Hz, Q = 2tt x 3kHz, k = 0.74. 
The frequency vo is the inverse period of the system, and Av 
represents the global phase rotation. From the data in this 
figure the values vo = 6.65kHz and Av = —1.78 kHz can be 
determined, giving the period as T — 150/xs. 



As an example of this process, we consider a two-level 
condensate with physical parameters corresponding to 
those used in the metastable helium BEC experiment 
by Dall et al. [35]. While a more detailed analysis of 
this experiment using our theory will be done in Sec- 
tion [Vl] we note that it fulfils the criteria we require — a 
multi-state, high aspect ratio condensate with the states 
possessing dissimilar scattering lengths and coupled by 
radio frequency radiation. Specifically we take a parti- 
cle number of N = 2 x 10 6 , radial and axial trapping 
frequencies of uj r — 2tt x 1020Hz and u z = 2tt x 55Hz re- 
spectively, choose a coupling field with a Rabi frequency 
of CI = 2tt x 3kHz, and nonlinearity mismatch parameter 
of k = 0.74. 

We simulate the evolution of the mean field of such a 



BEC using Eqs. ( 19 ), and obtain a time series of the two 



condensate fields &j(t) at the center of the condensate. 
We then Fourier transform these time series and obtain 
a power spectrum for the two fields, which is plotted in 
Figure [2] The periodic behaviour of the fields is clearly 
visible, and we can easily extract the period of T = 150/is 
and the frequency offset of Av = —1.78kHz. 

The period and energy offset determined, it remains 
to calculate the monodromy matrix A^(k) from which 
the Floquet exponents may be derived. This is achieved 
by numerically solving the related matrix problem to 
Eq. pil from t — to t = T, as can be seen from Eq. (|34|. 



Noting that the matrix H(k) only depends on k = |k| , the 
solutions for the Floquet exponents £(fc) = ui(k) + i-f(k) 
are shown in Figure [3] 

In the limit that the mean field is time-independent 
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(a degenerate case of periodicity) , the Floquet exponents 
£(k) are related to the eigenvalues A(k) of the matrix H 
by £(k) = \(k)/h. It was previously stated that eigenval- 
ues of 'H(k) with non-zero imaginary components would 
be unstable, this is also true for the Floquet exponents 

e(k). 

The temporal periodicity of the system implies that the 
real components (k) of the Floquet exponents are only 
uniquely defined modulo 2inso. This is because any eigen- 
value A of the monodromy matrix A^(k) corresponds to 
infinitely many Floquet exponents, 



A = cxp [— i (lu + 2tttivq 



■ ij) T] — exp [— i (u> + ij) T] . 

(36) 



This does not hinder our understanding of the system as 
our primary interest is in the stability of the condensate 
to excitations, which is determined by the 7i(k). 

The normalised eigenvectors of the system are not al- 
ways annihilation or creation operators; as is discussed 
in Appendix [Aj this is only true when the Floquet ex- 
ponents are purely real. When the Floquet exponents 
have a non-zero imaginary component, annihilation and 
creation operators can be constructed from linear combi- 
nations of the eigenvectors. Not being eigenvectors, these 
operators will therefore have non-trivial evolution. As is 
shown in Appendix [XJ Floquet exponents with non-zero 
imaginary parts come in pairs of the form w(k) ± «7(k). 
From the corresponding eigenvectors to these Floquet ex- 
ponents the bosonic annihilation operator A(k, t) can be 
formed, which evolves as 

A(k, nT) = e mu ^ T (sinh(n7(k)T)A't(-k, 0) 

+ cosh(n 7 (k)T)A(k,0)) , (37a) 

A'(k,nT) = 

+ cosh(n7(k)T)A'(k,0)) , (37b) 



e -»M k ) T (sinh(n 7 (k)r)A t (-k, 0) 



where n is a positive integer, and where A'(k, t) will be 
equal to A(k, t) in some circumstances as discussed in 
Appendix [A] Due to the exponential growth in Eq. (371, 
the mode corresponding to A(k) represents a dynamical 
instability of the condensate. 

The behaviour of the dynamical instabilities is gov- 
erned by Eq. (37) only while the unstable modes have a 



small occupation compared to the condensate, and scat- 
tering between the unstable modes can be neglected. Fig- 
ure [3^c) shows the results of a truncated Wigner simu- 
lation of the Hamiltonian Eq. ( |15[ ), which is in excellent 
agreement with the location of the dynamical instabil- 
ities as determined by the Floquet exponents [see Fig- 
ure [3^b)]. For later times, there is an additional mode 
undergoing growth, k « 7.5 x 10 5 m _1 . This mode is the 
result of scattering between the dynamical instabilities, 
a process neglected the perturbative approach taken in 
this section. The depletion of the condensate mode has 
also been neglected in this analysis. 
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FIG. 3: (color online) Illustration of the Floquet exponents for 
H(k) and comparison to a corresponding truncated Wigner 
simulation for fl = 2tt x 3kHz. Upper figure (a) displays the 
real part uui of the Floquet exponents. As the uji can only be 
determined up to a multiple of (see Eq. ( 36 )), they have 
been reduced modulo 2-kv§ into the range \--kvq, 7t^o]- The 
middle figure (b) shows the imaginary part 7; of the Floquet 
exponents, which indicate an instability for the corresponding 
wave number when they are non-zero. Note that what appears 
to be a horizontal line at 7^ = is not an axis, but a plot of 
the 7i(fc), which are mostly zero. Thus at k = 8.26 x 10 5 m _1 
all four modes are unstable, but at k = 2.25 x 10 6 m _1 only 
two modes are unstable. As discussed in the Section |VIII| 
the two-mode regime is more suitable for experimental detec- 
tion of entanglement due to the entanglement being between 
fewer modes. Lower figure (c) shows the results of a ID trun- 
cated Wigner simulation corresponding to the system under 
consideration in (a) and (b). The truncated Wigner simula- 
tion exhibits growth in the same modes predicted from the 
results of the perturbative analysis shown in (b). The trun- 
cated Wigner results shown are the average of 500 realisations. 
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In summary, the procedure used to find the Floquet 
exponents of "H(k), and hence the stability of the con- 
densate to excitations is: 



1. Numerically solve Eq. (19) with fi — for a long 



time t ^> T where T is the periodicity of the solu- 
tion. 

2. Perform the temporal Fourier transform of the so- 
lutions obtained for ^i(t) to accurately determine 
the period T and the value of the energy offset jU re- 
quired to cancel any global phase evolution to make 
the wave functions themselves periodic. 



3. 



Using the calculated period and energy offset, nu- 
merically solve the related matrix problem given 
by Eq. ( 24 ) for a range of values of k to obtain the 
monodromy matrix A4(k). In this calculation the 
matrix A(t) in Eq. Q is -iH(k,t)/h. 



4. Calculate the eigenvalues of A4(k) and 
determine the Floquet exponents £j using 
Xi = cxp(— i&T). The real parts of these Flo- 
quet exponents give the energy spectrum, with 
non-zero imaginary parts giving the growth rate 
for the corresponding instability. 



VI. AGREEMENT OF THEORY WITH 
EXPERIMENT 

We have shown that if we apply a coupling between the 
two levels of a two-state condensate, then, provided the 



nonlinearity mismatch parameter k given by Eq. ( 18 ) is 
not equal to one, there exist instabilities which can result 
in excitations that grow exponentially at well-defined mo- 
mentum values fik. These excitations, however, will only 
become amplified provided their de Broglie wavelength is 
smaller than the extent of the condensate in the direction 
of k — if the entire condensate covers only a fraction of 
a wavelength, clearly no amplification can take place. 

To demonstrate the validity of our theory, we apply it 
to the experiment of Dall et al. [32] . In this experiment 
metastable helium atoms in the mp — +1 Zeeman state 
were held in magnetic trap and evaporatively cooled to 
form a BEC. An RF coupling field was then applied, flip- 
ping the atoms from the mp = +1 magnetically trapped 
state to the mp = untrapped state, forming an atom 
laser. These untrapped atoms fell from the trap onto a 
position sensitive detector 4cm below and the resulting 
two-dimensional beam cross section was imaged. In addi- 
tion to the expected atom laser profile, they also observed 
two anomalous beams of atoms that had been ejected in 
the weak trapping direction of the condensate. 

One complication is that metastable helium in the 
F = 1 manifold possesses three possible Zeeman states, 
mp = +1, 0, —1, and our theory was developed for two- 
state systems. As the mp = — 1 state is antitrapped, 
any atoms in this state leave the condensate very rapidly 



due to the combined effects of both the mean field re- 
pulsion and the magnetic field gradient. For the ex- 
periment in question, a classical particle in the centre 
of the condensate under the influence of the same effec- 
tive potential experienced by the mp = — 1 atoms would 
reach a momentum equal to the momentum width of the 
condensate (and hence no longer be able to couple to 
the stationary atoms in the middle of the condensate) in 
~ 80/^s, significantly shorter than the inverse Rabi fre- 
quency of ~ 300/Lts. This indicates that the population of 
the mp = — 1 state will be low, and have minimal influ- 
ence on the dynamics of the system, and consequently we 
make an approximation where we ignore the mp = — 1 
state. 

The relevant experimental parameters were a particle 
number of N = 2 x 10 6 , and radial and axial trapping 
frequencies of u> r = 2ir x 1020Hz and oj z = 2ir x 55Hz 
respectively, and the RF outcoupling field corresponding 
to maximum anomalous particle production had a Rabi 
frequency of fl — 2ir x 3kHz. Metastable helium has 
s-wave scattering lengths of ai t \ — a^o = 7.51nm and 
a o,o = 5.56nm |54j . where a,i.j is the scattering length 
between the two Zeeman states mp = giving a non- 
linear mismatch parameter of k — 0.74. 

We took two approaches to numerically modelling this 
system. The first was to simply solve the coupled Gross- 
Pit aevskii (GP) equations, which gives the mean field 
behaviour of the BEC and outcoupled atoms. The GP- 
equations, however, are a mean field approximation, and 
our theory involves physics beyond the mean field. Con- 
sequently the second approach was to numerically simu- 
late the system using a phase-space method, specifically 
the truncated Wigner approximation (TWA). 

In both cases we simulated the full system in three 
dimensions in a fully multi-modal way, thus including 
spatial density variations across the condensate. This 
contrasts with the simpler numerical calculation carried 
out in Section [Vj where we assumed a density that was 
constant in space but not time. As any real condensate 
has spatial density dependence, this is an important test 
of the validity of our model. 

The results of these simulations are plotted in Figure |1J 
showing both the momentum distribution of the outcou- 
pled atoms as well as the cross-sectional density profile 
of the atom laser on the detector. It is clear that the GP 
simulation shows no significant momentum components 
in the axial (weakly trapping, long condensate axis) di- 
rection. The TWA simulation on the other hand, shows 
significant specific momentum components have been ex- 
cited along the axial direction. 

We have also plotted where we expect these momen- 
tum components to lie based on the perturbation analysis 
in Section [V] The dynamical instability with the largest 
growth rate illustrated in Figure [3] at k = 2.25 x 10 6 m _1 
is in good agreement with the highest growth rate insta- 
bility in Figure [4^c) observed at (k r rj 0, k z « ±2.4 x 
10 6 m -1 ), and the position at which the particles corre- 
sponding to this instability land on the detector is in good 
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FIG. 4: (color online) Simulation results for outcoupling from 
the centre of the condensate. Left figures (a) and (c) plot the 
momentum density of the untrapped state at t = 2 ms as 
a function of the radial (k r ) and axial (k z ) wave numbers. 
Right figures (b) and (d) plot the normalised density profiles 
that would be observed on the MCP detector 4 cm below 
the condensate due to outcoupling from the condensate for 
t = 6 ms. The upper figures (a) and (b) plot the results of 
a GP model where the mp — —1 state is assumed negligible, 
lower figures (e) and (f) correspond to a Truncated Wigner 
simulation for the same system averaged over idealisations = 4. 
The wave numbers of the three fastest-growing instabilities 
predicted by the perturbation analysis performed in Section|V] 
are marked in (c). The weak axis is in the vertical direction 
in all figures. 



agreement with the images from the experiment [52] . 

The absence of the effect in the GP simulations and the 
presence of the effect in the TWA simulations indicates 
that dynamic instabilities causing amplified excitations 
are indeed present in such systems, and can be experi- 
mentally created and measured. 



VII. EXISTENCE OF ENTANGLEMENT 

While the production of oppositely directed matter 
wave beams though dynamic instabilities in a condensate 
is interesting, it is also worth noting that these beams 
should exhibit EPR entanglement on formation. 

EPR entanglement was initially proposed from a first- 
quantized point of view |22j . Two individual particles 
are created with entanglement between two conjugate ob- 
servables — position and momentum, for example, or two 
orthogonal spins — and then the particles are separated. 
Due to the entanglement, the measurement of one observ- 
able on one particle forces the other to instantaneously 
adopt a predictable value for the same observable, no 



matter how large the distance between the two particles. 

In modern implementations, rather than considering 
discrete particles, one usually considers continuous vari- 
able entanglement between the quadratures of atomic or 
optical fields instead. Provided the correlation functions 
between the quadratures violate a specific inequality, the 
fields are EPR-entangled [26] . 

The standard procedure to create such fields in the 
optical regime is to use parametric down conversion [55] . 
Here a single photon of frequency 2w interacts with a non- 
linear crystal and is split into two photons of frequency u> 
in the degenerate case and two photons with frequencies 
lui +UJ2 = 2w in the nondegenerate case. The two gener- 
ated photons are EPR correlated, and it is this process 
that has been the workhorse of almost all optical EPR 
tests. 

As discussed in Section [Vj our scheme generates quasi- 
particles at specific momenta, and it is possible to derive 
the equations of motion for the bosonic operators de- 
scribing these quasiparticles. Restated, these equations 
of motion are 

A(k,rcT) = e m " (k)T (sinh(n 7 (k)T)A' t (-k, 0) 

+ cosh(n 7 (k)r)A(k,0)) , (38a) 
A'(k,nT) = e - inu(k)T (sinh(n7(k)T)At(-k,0) 

+ cosh(n 7 (k)T)A'(k,0)) . (38b) 



These equations describe the exponential amplification 
of two distinct modes, with each mode existing as a su- 
perposition of the |1) and |0) states. Equations (38) are 



identical to those describing optical two-mode paramet- 
ric down-conversion [55], with the Floquet exponent 7 
playing the role of the strength of the crystal nonlinear- 
ity. Just as these equations lead to EPR entanglement 
for photons, here they lead to EPR entanglement be- 
tween the quasiparticle excitations. We note that Equa- 
tions ( 38 ) do not strictly give the continuous time evolu- 
tion of the atomic fields, but rather are only stated for 
times that are multiples of the Floquet period T. How- 
ever as the Floquet exponents govern the long term be- 
haviour of the system, the entanglement of the atomic 
system will asympotically approach that of the analogous 
optical system, even for times t =/= T. 

The entanglement produced by the evolution described 
by Eq. ( 38 ) is between the oppositely-propagating modes 
A(k) and A'(— k). Strictly, Eq. (38) only applies when all 



the real parts of all four Floquet exponents are of equal 
magnitude, and the imaginary components are of equal 
magnitude as is the case near k — 8.26 x 10 5 m" 1 in Fig- 
ure [3j When only two Floquet exponents have non-zero 
imaginary components (or the real parts are of different 
magnitude), A'(k) = A(k), and the entanglement will be 
between A(k) and A(— k). Equations (38) are consistent 
in this case as either w(k) = or w(k) = ttisq, giving 
exp (inw(k)T) = ±1. 
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As the A and A' are superpositions of the |1) and |0) 
states, the entanglement will only be measured in the ap- 
propriate superpositions of these states. For the case in 
which only two Floquet exponents have non-zero imag- 
inary components (i.e. A(k) = A'(k)) there will only 
be one superposition of the |1) and |0) states that will 
be formed at that momentum (the other being stable 
will not undergo exponential growth), and it will be the 
oppositely-propagating mode of this superposition with 
which it will be entangled. In the case in which all 
four Floquet exponents have equal (non-zero) magnitude 
imaginary components and equal magnitude real com- 
ponents, the entanglement will be between orthogonal, 
oppositely- propagating superpositions of the |1) and |0) 
states. For example, if A(k, 0) were purely the |1) state, 
it would be entangled with the |0) state propagating in 
the opposite direction. 

It should also be noted that this analysis holds true 
both for the case where the A and A' are spontaneously 
seeded, and for the more general case in which the A and 
A' are seeded by different coherent states. 



VIII. DISCUSSION AND EXPERIMENTAL 
PRACTICALITY 

We have demonstrated that our scheme can create par- 
ticle beams along the long axis of a BEC, and shown that 
the generation of these these beams has the same mathe- 
matical form as optical non-degenerate parametric down 
conversion, resulting in EPR-entangled pairs of excita- 
tions with oppositely-directed momenta. In addition, an 
experiment satisfying the criteria of our scheme has been 
shown to exhibit cones of particles ejected along the long 
axis of a BEC, with momenta generally agreeing with our 
predictions. 

It is clear that since the entangled modes are excita- 
tions described by the A(k) operators, it is not a simple 
case of, say, an mp = atom being entangled with an 
mp = 1 atom with opposite momentum, as the entan- 
gled A(k) modes are each superpositions of mp — 1 and 
mp = states. This is a problem in an experiment where 
the two states are not confined in the same potential, as 
in the metastable helium experiment discussed earlier. In 
such a case, the mp = 1 atoms remain magnetically con- 
fined, while the mp = atoms can leave the condensate. 
However, if the mp = 1 components of the entangled 
modes are outcoupled faster than the time to undergo 
half an oscillation in the trap and reverse their momenta 
(~ 9ms in the case of the helium experiment), then it 
should be possible to at least observe number difference 
squeezing between the atoms emitted in opposite direc- 
tions. Even if the mp — 1 atoms have enough momen- 
tum to escape, however, the fact that they see a different 
potential to the mp = atoms leads to considerable dif- 
ference in their dynamics. 

One solution would be to use optically trapped atoms, 
with the Zeeman substates split via an external magnetic 



field. This will ensure that each substate will see the same 
confining potential. The difficulty with this approach is 
most of the conveniently accessible systems have at least 
three Zeeman states, ranging over {—j, —j + 1, . . . j}. If 
RF radiation is used to couple these states, then since 
the energy splitting between the states is (usually) lin- 
ear, all three states will become populated, breaking our 
two-level assumption. One might think that in such 
a three state system one could couple mp = —1 and 
mp = +1 via a Raman transition, but such two-photon 
Amp — 2 transitions within the same F manifold are 
heavily suppressed [8]. A better approach would be to 
use magnetically trapped BECs in two separate hyper- 
fine manifolds, and couple two specific Zeeman states — 
\F = l,mp = —1) and |F = 2,m F — +1), for example. 
Another approach could be to use the quadratic Zeeman 
shift to remove the degeneracy of the energy splittings 
in an F manifold. Using a sufficiently large bias field in 
conjunction with an optical trap would make RF cou- 
pling between the m = 1 and m = levels of an F = 1 
manifold effectively a two-level system. 

If a two-level system cannot be obtained, a three-level 
system could still be used. This is not a serious prob- 
lem provided all atoms begin in a single trapped state 
(whichever one of the mp — ±1 states is trapped in the 
magnetically trapped case, or either of these states in the 
optical trapping case), couple to the mp = state, and 
ensure the coupling is weak enough that the atoms es- 
cape more quickly than the period associated with the 
Rabi frequency of the coupling. This will ensure that al- 
most all of the population remains in two levels, and the 
system behaves like a two-level system. As discussed in 
Section |VI[ this approximation describes the case for a 
magnetically trapped metastable helium BEC quite well, 
despite its three levels. 

Although the particles forming the entangled beam will 
usually have enough energy to escape the trapping poten- 
tial, it is possible that for some choice of atomic species, 
atomic density, and trapping geometry the momentum 
at which the effect occurs may not be be high enough. 
In this case the experiment could be run in pulsed mode, 
dropping the trapping potential after the entangled ex- 
citations have formed, and allowing them to escape. As 
these entangled particles leave along the weakly trap- 
ping axis, they will be well-resolved from the bulk of the 
condensate which expands most rapidly along the tight 
trapping axis. 

The actual measurement of EPR quadrature entangle- 
ment requires the existence of a local oscillator to turn 
the quadrature entanglement into number correlations. 
Local oscillators are difficult to obtain for atomic fields; 
generally the only option is to use the source condensate 
itself. In this case the condensate would be split in two, 
with one of the condensates used to generate the entan- 
gled beams and the other used to generate the atomic 
beams that will form the local oscillator. As the conden- 
sates will initially only be displaced by a few hundred 
microns, the beams will overlap well in the far field al- 
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lowing reasonable mode matching. 

The most experimentally promising regime in which 
to observe entanglement would be to consider atomic 
beams for which only two Floquet exponents have non- 
zero imaginary components (this is observed, for exam- 
ple, near k — 2.25 x 10 6 m" 1 in the He* experiment). 
In this regime the two counter-propagating local oscil- 
lators will have the same superposition of |1) and |0) 
states, presenting less of a challenge than producing 
oppositely-propagating atomic beams in orthogonal su- 
perpositions. Additionally the correct superposition of 
|1) and |0) states can be determined experimentally up 
to a relative phase difference by measuring the fraction 
of |1) and |0) atoms in the produced atomic beams. If 
the correct superposition were not known in advance it 
would otherwise be necessary to scan over all of them, 
which is a two-dimensional space. Knowing the ratio of 
|1) and |0) atoms would reduce the search space to a 
single dimension. 

The regime in which all four Floquet exponents have 
non-zero imaginary components of equal magnitude (and 
equal magnitude real components) is interesting as two 
pairs of entangled counter-propagating beams are pro- 
duced. Measuring this entanglement is however a greater 
challenge as the local oscillators are more complicated, 
and it would require searching over a larger space to find 
the entanglement unless highly accurate theoretical pre- 
dictions for the experiment were available. 

The generation of the beams also relies on the correct 
choice of atomic species. The effect depends on a non- 
linearity mismatch, requiring that the s-wave scattering 
lengths in the system are not all identical, so that k^I. 
This is not necessarily easy to achieve. 87 Rb, for exam- 
ple, is a popular and well-understood atom in the BEC 
community, but if the F = 1 hyperfine manifold is used 
for both Zeeman sublevels to try and remain close to 
two-state behaviour, then k = 1.002 56], and even if we 
consider using one level from the F = 2 manifold and 
one in the F = 1 manifold and coupling the two via an 
optical Raman transition, the phase mismatch is at best 
a few percent [57] . 



IX. CONCLUSION 

We have investigated the evolution of a two-state BEC 
where the two states |0) and |1) are coupled in some 
way, and the s-wave scattering lengths between (|0), |0)), 
(|1), |1)) and (|0), |1)) are not all identical. Our analy- 
sis demonstrates that the mean field background of such 
a system exhibits periodicity, in that the same physical 
state occurs repeatedly with a fixed period. 

To investigate the stability and energy excitation spec- 
trum of the system, we performed a perturbation analy- 
sis beyond the mean field. As the mean field background 
about which we are considering perturbations is time- 
dependent, the standard approach would be to apply the 
Bogoliubov-de Gennes theory, which would result in a 



purely numerical solution to the dynamics. As our sys- 
tem is time-dependent but periodic, however, we are able 
to apply Floquet 's theorem to obtain analytic insights 
into the system. 

From the period of the mean field evolution, we can 
extract the monodromy matrix describing the dynamics, 
the eigenvalues of which allow us to calculate the Flo- 
quet exponents for the system. These Floquet exponents 
are generally complex, indicating dynamic instabilities 
within the system. Specifically, the real part of the expo- 
nents describe the energy spectrum of excitations within 
the condensate, while the imaginary part corresponds to 
exponential growth of the excitations. As the exponents 
are momentum-dependent, this results in the amplifica- 
tion of excitations at very specific momenta. 

The evolution of the annihilation and creation opera- 
tors of the excitations are identical to those correspond- 
ing to the annihilation and creation operators of the pho- 
tons in optical degenerate parametric down-conversion. 
Optical parametric down-conversion is a well-studied and 
often-employed scheme to generate EPR-entangled pairs 
of photons and, just as in the optical case, our scheme 
results in EPR entanglement between two particle excita- 
tions with the same momentum magnitude but opposite 
direction. 

As the amplification can only occur if there exists sub- 
stantial condensate population within a few de Broglie 
wavelengths of the initial excitation, it is possible to cre- 
ate a situation when the entangled excitations are emit- 
ted directionally. This can be done by employing a high 
aspect ratio condensate, where the Thomas-Fermi radius 
in the tightly trapped axis is comparable to or less than 
the de Broglie wavelength of the lowest momentum exci- 
tation, resulting in exponential gain predominantly along 
the long axis of the condensate. The narrower the con- 
densate is, the more directional the entangled beams can 
become. 

The system we have analyzed in this paper is not diffi- 
cult to arrange experimentally, and we have shown that 
an experiment matching our criteria resulted in particle 
emission along the long axis of the condensate at momen- 
tum that matches our predictions |32j . That experiment 
was not optimized for the production of entangled parti- 
cle beams, however, and, as we have discussed, a better 
scheme would involve an optically rather than magneti- 
cally trapped condensate. 

We would like to thank Joseph Hope and Karen 
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for Quantum Atom Optics and ARC Discovery Project 
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Appendix A: Evolution of the excitations 

The theory of elementary excitations in unstable Bose- 
Einstein condensates has been considered before |42) . 
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In this appendix, restrictions on the Floquet exponents 
for the system under consideration are obtained from 
which the equations of motion for the instabilities may 
be solved. This is achieved by applying the methods 
discussed in [52] and consideration of the symmetries of 
the system Hamiltonian. The excitation dynamics is de- 
scribed by a temporally periodic but spatially homoge- 
nous and isotropic Hamiltonian. This Hamiltonian can 
be obtained by making a perturbative expansion of the 
Hamiltonian Eq. ( 15 ) considered in Section III about the 



time-dependent (but periodic) mean-field. 

Consider a general quadratic, spatially-homogeneous 
Hamiltonian H(t) of period T in terms of the opera- 
tors <f>i(x,t) which obey the usual equal-time bosonic 
commutation relations. Due to conservation of mo- 
mentum, the Hamiltonian H(t) when written in the 
Fourier basis can only contain terms of the form 
$(k,tU.(M), #(k,t)#(-k,t) and 4 (k,<U.(-k,t), 



while terms of the form </>J(k, t)(j>j{— k, t), <^>J(k, *)0}(k, t), 
and ^(k, t)(f>j(k, t) cannot occur. As the Hamiltonian 
H (t) is homogeneous by assumption, the operator equa- 
tions of motion will take their simplest form in a Fourier 
basis. In this basis the equations of motion for the oper- 
ators 0j(k, t) can be written in matrix form 

*ft^T(k,t)=«(M)T(M), (Ala) 



dt 



T(M) = (^(k,i) #(-M) &(M) ••■ 

(Alb) 



where the matrix "H(k,t) obeys 

W(k,t + T) =H(k,t), 
H(k,t)=H(-k,t), 



(A2) 
(A3) 



where the last equality holds because H (t) is isotropic. 

If the Hamiltonian H(t) were not time-dependent, 
H(k, t) could be diagonalized to find the (potentially 
complex) eigenvalues flj (k) and corresponding operators 
Q j (k, t) which would evolve as 



(A4) 



where the Q j (k, t) need not obey boson commutation re- 
lations. The real parts of the eigenvalues of H(k, t) would 
give the excitation spectrum of the Hamiltonian H{t) 
with non-zero imaginary components giving the growth 
rates of the corresponding unstable mode. 

In the case of a periodic matrix W(k, t) it is instead the 
monodromy matrix A4(k) (see Section M) that we wish 
to diagonalize. The monodromy matrix M(k) satisfies 



T(k,nT) =M(k)"T(k,0), 



(A5) 



where n is a positive integer. In place of Eq. (A4) we 
seek the operators Qj (k, t) that obey 



Q 3 -(k > r)=A J -(k)Q J -(k,0) > 



(A6) 



where the Qj(k, t) are denned by 

Qj(k,t) = c\(k)T(k,t), 



(A7) 



for some vectors Cj(k), where cKk) denotes the conjugate 
transpose. 



Using definitions Eq. ( A5 |-Eq. ( A7), it follows that the 



Xj (k) and (k) are respectively the eigenvalues and left 
eigenvectors of M(k), 

Oj(k,T) = cj(k)f (k,T) = cj(k)7W(k)T (k,0) 



Qj (kj T) = A,(k)4(k,0) - A 3 -(k)ctt(k,0) 
=> c){k)M{k)± (k,0) - A 3 -(k)ctt(k,0) 
=> c](k)X(k) = A,(k)cJ(k), 



(A8) 
(A9) 
(A10) 
(All) 



where Eq. JAlTJ follows from Eq. ( |A10[ ) as the compo- 
nents of T(k,0) are linearly independent operators. 

The operators Qj (k, t) are not necessarily bosonic an- 
nihilation or creation operators. To determine the condi- 
tions under which they are, we consider their Hermitian 
conjugates Qj(k,t). As every operator in Y(— k,i) is 
the Hermitian conjugate of an operator in T(k, t), the 
Qt (k, t) can be written as 



Q }(k,t) = d (k)Y(-k,i) 



(A12) 



Tor some vectors dj(k). It follows from Eq. (A6| that the 
Q\(k,t) will obey 



Qt(k,T) = A*(k)Qj(k,0). 



(A13) 



The commutators of the Qj (k,t) will be constant as 

they are constant linear combinations of the 0^(±k, t), 
the commutators of which are themselves constant. Us- 
ing this requirement gives 

Q 2 (k,T), g](k,T)] = [Ai(k)O i (k,0) ) A*(k)Qt(k,0) 

(A14) 

= Ai(k)A*(k) [^(k.0), 0J(k,0)" . 

(A15) 



For Eq. (A15) to be true either A 4 (k)A*(k) = 1 or the 

two operators commute. Specifically, for Qi(k, t) to be 
an annihilation or creation operator it is required that 
A*(k) _1 = Aj(k). In terms of the Floquet exponents (see 

Section V) £j(k) = — In A^(k), this requirement becomes 

£i(k) = Wj(k). Hence it is only for purely real Floquet 
exponents that the eigenvalues of Ai(k) correspond to 
bosonic annihilation or creation operators. Note that in 
the degenerate case in which T-L(k, t) is time-independent, 
the Floquet exponents £(k) are equal to the eigenvalues 
f2(k) of H(k,t). Hence the real components of the Flo- 
quet exponents are related to the excitation spectrum 
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and non-zero imaginary components are related to the 
existence of instabilities. 

Generally the Floquet exponents £j may have a non- 
zero imaginary component. In this case the Qi(k, t) will 
not be bosonic annihilation or creation operators, al- 
though such operators can be constructed from linear 
combinations of the Qj(±k, i). Before constructing such 
operators, we first consider the restrictions on the possi- 
ble eigenvalues of M (k) . 

First it is noted th at if A, (k) is a n eigenvalue of 
M{k), then from Eq. jAl2| -Eq. fAl3| A*(k) must be 
an eigenvalue of M.(— k). However, from the reflection 



symmetry of H(k) defined by Eq. (A3 1 we have that 
M(— k) = M(k) and hence A*(k) must also be an eigen- 
value of M (k) . 

Secondly, not all of the operators Q^\k,t) can com- 
mute. As the Q^p (k, t) form a complete basis over the 
same space as the $p (k, t) which themselves do not all 
commute, for every operator Qi (k, t) there must be at 
least one other operator with which it does not commute. 
From Eq. ( A15 ) then follows the requirement that if Aj(k) 
is an eigenvalue, A* (k) must also be an eigenvalue. 

Combining these two requirements gives a consistency 
condition for the eigenvalues of M.{k): if A is an eigen- 
value, A*, A -1 , and A* -1 must all be eigenvalues. These 
conditions can be met using 1, 2 or 4 distinct eigenvalues 
ofM(k). 

For the degenerate case in which all of A, A*, A -1 , 
and A* -1 are equal, the eigenvalue A = ±1. The corre- 
sponding Floquet exponent is £ = or £ = ttvq where 
Uq = T^ 1 . This is not an interesting case and does not 
occur in Figure [3j 

There are two ways that two distinct eigenvalues can 
be used to satisfy the consistency condition for the eigen- 
values. The first possibility is that A = A* -1 (with A* 
being the second eigenvalue). In this case the Floquet 
exponents are £ = ±w. In this case, the operators corre- 
sponding to the eigenvalues are bosonic annihilation or 
creation operators as shown above. The second possi- 
bility is that A = A* (with A -1 being the second eigen- 



value). In this case the Floquet exponents are £ = ±ij 
or £ = 7r^o±*7- This situation is seen in Figure [3] around 
k « 1.5 X lC^m" 1 and k 2.25 x 10 6 m -1 . 

The final possibility is that four distinct eigenvalues 
are used to satisfy the consistency condition. In this case 
the four eigenvalues A, A*, A -1 , and A* -1 are different 
and the corresponding Floquet exponents are £ = u> ± 17 
and £' = —lj ± i"f. It is this situation that occurs in 
Figure |3] around k w 0.75 x 10 6 m -1 . 

In summary, the Floquet exponents with non-zero real 
parts come in pairs £(k) = u>(k) ± i"f(k). From the oper- 
ators corresponding to these pairs of exponents bosonic 
annihilation and creation operators can be constructed. 

Consider the eigenvalues A = e r+l< ^ and A' = e~ r+l< ^, 
and the corresponding operators Q(k,t) and Q'(k,t). 
From these operators we define the following two op- 
erators which will respectively be shown to be bosonic 
annihilation and creation operators, 



(A16) 
(A17) 



A 1 (k,t) = -^=[Q(k,t) + Q'(k,t) 
A 2 (k,i) = -L(Q(k,i)-g'(k,t) 



As A A'* = 1, Q(k,t) and Q't(k,t) will not commute. By 
appropriate rescaling of the operators, we can define their 
commutator to be 



Q(k,t), Q'\k,t)] =1 



(A18) 



This choice defines the value of the other non-zero com- 
mutator, 



Q'(k,t), Qt( k)t )l = i 



(A19) 



From these two commutators it can then be shown that 
Ai(k, i) obeys the commutation relations appropriate for 
an annihilation operator, while A2(k, t) obeys the com- 
mutation relations for a creation operator. For example, 



AitM), A\(k,t) = - Q(k,t) + g / (k,i),Qt(k,t) + Q't(k,t) 



1. 



(A20) 



Defining A^(— k, t) = A|(k, i), the evolution of the operators Ai(k, t) and A^(k, t) can now be determined. Ai(k, t) 
evolves as 



Ax(k,nT) = (Q(k,nT)+Q'{k,nT) 

3 nr+m ^Q(k,0) + e - nr+, ' : " 



V2 
1 

71 



Q'(k,0)) 



Jncji 



V2 



( e «r _ e -nr) Q( k , 0) - Q'(k, 0) + - (e nr + e~ nr ) -= Q(k, 0) + Q'(k, 0) 



V2 



(sinh(n7T)Ai t (-k, 0) + cosh(?i7T)Ai(k, 0)) , 



(A21) 
(A22) 

(A23) 
(A24) 
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where n is a positive integer. Similarly, A'^k, t) can be 
shown to evolve as 

A;(k,nT) = e- m " T ^sinh(n 7 r)A t 1 (-k,0) 

+ cosh(n 7 T) A; (k, 0) ) . (A25) 



Finally, if the Q(k, t) and Q'(k, t) operators correspond 
to only two distinct eigenvalues (i.e. A = A*), then by the 
uniqueness of the eigenvectors A' ( k, t) = Ai(k, t) . For 
this case e mujT = ±1, making Eq. ( |A24| and Eq. ( |A25[ ) 
consistent. 



Appendix B: Proof of periodicity of the nonlinear 
Bloch equations 

In this Appendix we give a proof that the solutions to 
the nonlinear optical Bloch equations considered in Sec- 



tion III which describe the dynamics of the mean field of 
a two-level homogeneous condensate in the case that the 
scattering lengths for the two levels are not equal, are pe- 
riodic. We restate the nonlinear optical Bloch equations 



(201 here for convenience: 



d Pw = -if (1 - w)pio + i^w, 



dt 



d_ 
dt 



w = -4film{/9i }, 



(Bla) 
(Bib) 



where pio is the off-diagonal element of the density ma- 
trix, and w = pn — poo- 

These equations describe motion on the surface of 
a sphere, known as the Bloch sphere (see Figure [l}, 
with each point on the sphere corresponding to a dif- 
ferent physical state. A generalisation of the Poincare- 
Bcndixson theorem that applies to compact, connected, 
two-dimensional orientable manifolds (such as the sur- 
face of the sphere) states that every trajectory either 
approaches one or more fixed points, or approaches a 
periodic orbit |48j . The possibility of a space-filling tra- 
jectory is precluded as the surface of a sphere is not home- 
omorphic to a torus (see [IB]). Periodic trajectories and 
fixed points trivially approach themselves. 

We next show that there are no limit cycles in this 
system, i.e. there exist no trajectories that approach pe- 
riodic orbits which are not themselves periodic. 

We first write the equations of motion for the system 
in the usual spherical polar coordinates as 



de 

~di. 



2D sin ( 



— =2Vt cos (f> cot 9 — 2-9(1 — cos ( 



(B2a) 
(B2b) 



where we note that these equations are of the form 



r = — r x VE, 

h 



(B3) 



where E is the energy per particle previously defined in 
Eq. (21), and a hat is used to denote unit vectors. 



Eq. ( B3 1 shows that a given point moves in the direc- 



tion perpendicular to the gradient of E at a rate propor- 
tional to the magnitude of the gradient, and hence the 
energy E is conserved along any given trajectory. 

Assuming a limit cycle exists, there will exist a trajec- 
tory approaching the limit cycle which is not the limit 
cycle. For each point in the limit cycle x there exists an 
infinite sequence of points that approach x which are 
the intersection of a curve which intersects no trajectory 
tangentially (a transversal, see [15]) and the trajectory 
approaching the limit cycle. As the energy E is con- 
tinuous, we have lim E(si) — E(x), and therefore the 

i— >oo 

energy of the approaching trajectory must be the same 
as the energy of the limit cycle. Further, as energy is 
conserved along the trajectory, and all the points lie 
on a trajectory, E(si) = E(x.) for all points s^. Thus 
the derivative of E in the direction of the transversal is 

zero as lim — p— ^ = 0. It is already known that 

i->oo |s.; — x| 

the derivative of the energy is zero parallel to the limit 
cycle as energy is conserved along any trajectory. As the 
direction of the transversal and the limit cycle at x are 
linearly independent (a transversal intersects no trajec- 
tory tangentially) and the manifold is two-dimensional, 
the gradient of E at x must necessarily be zero. From 
(B3) this implies that x is a fixed point, which is a con- 



tradiction with x being part of a limit cycle. Our original 
assumption that a limit cycle exists is therefore false. 

It now remains to show that no trajectories approach 
fixed points (with the exception of the fixed points them- 
selves). As the possibility of limit cycles has been 
excluded, the generalised Poincare-Bendixson theorem 
shows that all remaining trajectories are periodic. 

To analyse the stability of the fixed points, it is first 
necessary to determine their location. It immediately fol- 
lows from (B2a) that fixed points will satisfy sin0 = 0, 
i.e. fixed points are restricted to the great circle in the 
x-z plane. To parameterise this circle by a single coor- 
dinate, we alter the usual ranges of the spherical polar 
coordinates such that <p £ [— § , § ) and 6 £ [— n, ir). With 
this definition cos 0=1 over the entire circle, instead of 
cos0 = ±1 over opposite halves. 

Using the tan-half substitution s 



tan 1(9, it can be 



shown from ( B2b I that all fixed points satisfy 



gs 3 = - s 4 ). 



(B4) 



The stability of (B2| about the fixed points may be 



used to demonstrate that there are no trajectories that 
approach the fixed points (with the exception of the fixed 
points themselves). Linearizing (B2) about the fixed 
points yields 



— 
~dt 

dt 



2tt5(t), 
~2^2 



(s 4 + 3)66, 



(B5) 
(B6) 
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where 68 and 5<p are the deviations from the fixed point 



and the identity (B4| has been used to simplify the re- 



sult. The eigenvalues of this linear system of differential 
equations are 



A = ± 



O 2 



(B7) 



which are pure imaginary as both Q, and s are real. This 
implies that there exist no trajectories near any fixed 
points of the system that asymptotically approach the 
fixed point. It may then be concluded that with the 
exception of the fixed points themselves, all trajectories 
of the system (Bl) are periodic, and Floquet's theorem 



may be applied in Section W\ 
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